
#Import the required libraries

from PIL import Image
import numpy as np
import pandas as pd
import os
import random
import warnings
warnings.filterwarnings("ignore")
import matplotlib.pyplot as plt
%matplotlib inline


#Datasets
vehicles_dataset = '/Users/prashantiyer/Documents/ML in Python/data/vehicles'
non_vehicles_dataset = '/Users/prashantiyer/Documents/ML in Python/data/non-vehicles'

# Create an empty list to store the 1D arrays
transformed_images = []
labels = []  # Create an empty list to store the labels

# Function to convert images and assign labels
def convert_images(dataset_path, label):
    images = []
    for filename in os.listdir(dataset_path):
        if filename.endswith(".png"):
            image_path = os.path.join(dataset_path, filename)
            with Image.open(image_path) as img:
                # Convert the image to grayscale
                grayscale_img = img.convert("L")  # "L" mode stands for grayscale
                grayscale_array = np.array(grayscale_img)
                transformed_image = grayscale_array.flatten()  # Flatten the 2D array to 1D
                images.append(transformed_image)
                labels.append(label)  # Append the label
    return images

# Load and flatten images from both datasets with labels
vehicles = convert_images(vehicles_dataset, label=1)  #label = 1; vehicles
non_vehicles = convert_images(non_vehicles_dataset, label=0) #label = 0; non-vehicles 

# Merge both the datasets
transformed_images.extend(vehicles)
transformed_images.extend(non_vehicles)

# Combine the images and labels into a single list and shuffle
combined_data = list(zip(transformed_images, labels))
random.shuffle(combined_data)


#Split the combined data back into images and labels
transformed_images, labels = zip(*combined_data)


# Convert to numpy array
grayscale_images = np.array(transformed_images)


#length of every image vector - 64x64
len(grayscale_images[1])

4096


#Min-Max Scaling on the image vector
from sklearn.preprocessing import MinMaxScaler

scaler = MinMaxScaler()

if len(grayscale_images.shape) == 1:
    grayscale_images = grayscale_images.reshape(-1, 1)

# Fit the scaler on your data and transform it
scaled_images = scaler.fit_transform(grayscale_images)

# If you want to retain it as a numpy array (optional)
grayscale_images = np.array(scaled_images)


#Count the number of vehicles and non-vehicle images in the dataset
from collections import Counter

Counter(labels).keys() # equals to list(set(words))
Counter(labels).values()

dict_values([8968, 8792])


#Function to reshape the image to their original pixel values and display

def plot_gallery(images, titles, h, w, num_images):
    """Helper function to plot a gallery of portraits"""
    plt.figure(figsize=(15,5))
    plt.subplots_adjust(bottom=0, left=.01, right=.99, top=.90, hspace=.35)
    for i in range(num_images):
        plt.subplot(1, num_images, i + 1)
        image = images[i].reshape((h, w))  # Assuming each image is hxw pixels
        label = "Vehicle" if labels[i] == 1 else "Non-Vehicle"
        plt.imshow(image, cmap='gray')
        plt.title(label)
        plt.axis('off')


#Display 10 images present in the dataset
plot_gallery(grayscale_images, labels, 64, 64,10)


%%time

import time
from sklearn.decomposition import PCA

#PCA for dimension reduction
start_time = time.time()

# Number of principal components (eigenvectors) to be considered
n_components = 50  

#Define the original dimensions of the image
h = 64
w = 64

# Create a PCA instance and fit it to your data
pca = PCA(n_components=n_components)
pca.fit(grayscale_images)

pca_data = pca.transform(grayscale_images)

reconstructed_data_pca = pca.inverse_transform(pca_data)

pca_time = time.time() - start_time
print(pca_time)

plot_gallery(reconstructed_data_pca, labels, 64, 64,10)

4.293466806411743
CPU times: user 29.1 s, sys: 981 ms, total: 30.1 s
Wall time: 4.47 s


def plot_explained_variance(pca):
    import plotly
    from plotly.graph_objs import Bar, Line
    from plotly.graph_objs import Scatter, Layout
    from plotly.graph_objs.scatter import Marker
    from plotly.graph_objs.layout import XAxis, YAxis
    plotly.offline.init_notebook_mode() # run at the start of every notebook
    
    explained_var = pca.explained_variance_ratio_
    cum_var_exp = np.cumsum(explained_var)
    
    plotly.offline.iplot({
        "data": [Bar(y=explained_var, name='individual explained variance'),
                 Scatter(y=cum_var_exp, name='cumulative explained variance')
            ],
        "layout": Layout(xaxis=XAxis(title='Principal components'), yaxis=YAxis(title='Explained variance ratio'))
    })

plot_explained_variance(pca)


%%time

import time
from sklearn.decomposition import PCA

start_time = time.time()

#Randomized PCA


# Define the dimensions of the original images (e.g., 64x64)
w = 64
h = 64

# Number of principal components (eigenfaces) to keep
n_components = 50  # You can adjust this number as needed

# Create a Randomized PCA instance and fit it to your data
rpca = PCA(n_components=n_components, svd_solver='randomized')
rpca.fit(grayscale_images)

# Transform the data into the PCA space
rpca_data = rpca.transform(grayscale_images)

# Reconstruct the images
reconstructed_data_rpca = rpca.inverse_transform(rpca_data)

rpca_time = time.time() - start_time
print(rpca_time)
plot_gallery(reconstructed_data_rpca, labels, 64, 64,10)

3.9093310832977295
CPU times: user 25.2 s, sys: 2.02 s, total: 27.3 s
Wall time: 3.95 s


plot_explained_variance(rpca)


from sklearn.metrics import mean_squared_error
import time

# Standard PCA
pca_error = mean_squared_error(grayscale_images, reconstructed_data_pca)
explained_variance_pca = np.sum(pca.explained_variance_ratio_)

# Randomized PCA
randomized_pca_error = mean_squared_error(grayscale_images, reconstructed_data_rpca)
explained_variance_rpca = np.sum(rpca.explained_variance_ratio_)


print(f"Standard PCA Explained Variance (for 50 components): {explained_variance_pca:.4f}")
print(f"Randomized PCA Explained Variance (for 50 components): {explained_variance_rpca:.4f}")
print(f"Standard PCA Reconstruction Error: {pca_error:.4f}")
print(f"Randomized PCA Reconstruction Error: {randomized_pca_error:.4f}")
print(f"Standard PCA Time Taken: {pca_time:.4f} seconds")
print(f"Randomized PCA Time Taken: {rpca_time:.4f} seconds")

Standard PCA Explained Variance (for 50 components): 0.8554
Randomized PCA Explained Variance (for 50 components): 0.8553
Standard PCA Reconstruction Error: 0.0065
Randomized PCA Reconstruction Error: 0.0065
Standard PCA Time Taken: 4.2935 seconds
Randomized PCA Time Taken: 3.9093 seconds


#Comparing PCA vs rPCA using Random Forest Classifier

from sklearn.model_selection import cross_val_score
from sklearn.ensemble import RandomForestClassifier  

classifier = RandomForestClassifier(n_estimators=100, random_state=42)  # You can choose an appropriate classifier for your problem
cv = 5

# Perform cross-validation for PCA
pca_scores = cross_val_score(classifier, pca.fit_transform(grayscale_images), labels, cv=cv, scoring='accuracy')
print("PCA Cross-Validation Scores:", pca_scores)
print("Mean Accuracy (PCA):", np.mean(pca_scores))

# Perform cross-validation for RPCA
rpca_scores = cross_val_score(classifier, rpca.fit_transform(grayscale_images), labels, cv=cv, scoring='accuracy')
print("RPCA Cross-Validation Scores:", rpca_scores)
print("Mean Accuracy (RPCA):", np.mean(rpca_scores))

PCA Cross-Validation Scores: [0.93524775 0.9321509  0.93412162 0.93778153 0.9420045 ]
Mean Accuracy (PCA): 0.9362612612612612
RPCA Cross-Validation Scores: [0.93721847 0.93130631 0.94031532 0.93975225 0.93834459]
Mean Accuracy (RPCA): 0.9373873873873875


#Store a single image to a variable to view the daisy process
image = grayscale_images[1].reshape((h, w))
image

array([[0.59215686, 0.59215686, 0.58039216, ..., 0.62745098, 0.63137255,
        0.61960784],
       [0.58823529, 0.59607843, 0.57254902, ..., 0.62745098, 0.63137255,
        0.63137255],
       [0.58431373, 0.58431373, 0.57647059, ..., 0.61960784, 0.61568627,
        0.62745098],
       ...,
       [0.47058824, 0.47843137, 0.4745098 , ..., 0.47843137, 0.4745098 ,
        0.47843137],
       [0.4745098 , 0.4745098 , 0.47843137, ..., 0.4745098 , 0.4745098 ,
        0.4745098 ],
       [0.47058824, 0.47843137, 0.47843137, ..., 0.4745098 , 0.4745098 ,
        0.47058824]])


from skimage.feature import daisy

# lets first visualize what the daisy descriptor looks like
features, img_desc = daisy(image, 
                           step=20, 
                           radius=20, 
                           rings=2, 
                           histograms=8, 
                           orientations=8, 
                           visualize=True)
plt.imshow(img_desc)
plt.grid(False)
plt.show()


# Function to tranform the features to daisy
def apply_daisy(row,shape):
    feat = daisy(row.reshape(shape), step=10, radius=20, 
                 rings=2, histograms=8, orientations=4, 
                 visualize=False)
    return feat.reshape((-1))

%time test_feature = apply_daisy(grayscale_images[1],(h,w))
print(test_feature.shape)
0.026 * len(grayscale_images)

CPU times: user 4.57 ms, sys: 859 µs, total: 5.43 ms
Wall time: 2.54 ms
(612,)

461.76


%time 
daisy_features = np.apply_along_axis(apply_daisy, 1, grayscale_images, (h,w))
print(daisy_features.shape)

CPU times: user 1e+03 ns, sys: 1 µs, total: 2 µs
Wall time: 2.62 µs
(17760, 612)


#Given the dataset is huge, we have proces the 40% of the sample in order to visualize and try finding a pattern.
#40% gives a good space for representing the a good range of features
import random
import seaborn as sns
from sklearn.metrics import pairwise_distances


sample_percentage = 0.40  # percentage of sub samples considered

# Calculate the sample size based on the percentage
sample_size = int(sample_percentage * daisy_features.shape[0])

# Randomly select sample indices without replacement
sample_indices = np.random.choice(daisy_features.shape[0], sample_size, replace=False)
sample_data = daisy_features[sample_indices]

# Calculate pairwise distances for the sample data
sample_distances = pairwise_distances(sample_data)

# Visualize the heatmap for the sample data
plt.figure(figsize=(10, 10))
sns.heatmap(sample_distances, cmap='viridis')
plt.title("Pairwise differences among sampled extracted features")
plt.show()


from sklearn.neighbors import KNeighborsClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score


# Split the data into training and testing sets
X_train, X_test, y_train, y_test = train_test_split(daisy_features, labels, test_size=0.2, random_state=42)

# Initialize a nearest neighbor classifier
knn_classifier = KNeighborsClassifier(n_neighbors=5)

# Train the classifier on the training data
knn_classifier.fit(X_train, y_train)

# Make predictions on the test data
y_pred = knn_classifier.predict(X_test)

# Evaluate the classifier's performance
accuracy = accuracy_score(y_test, y_pred)
print("KNN Accuracy:", accuracy)

KNN Accuracy: 0.9442567567567568


#Calculate the confusion matrix to know the cases that it does not classify correctly.

from sklearn.metrics import confusion_matrix
import seaborn as sns
import matplotlib.pyplot as plt

# Compute the confusion matrix
confusion = confusion_matrix(y_pred, y_test)

# Plot the confusion matrix as a heatmap
plt.figure(figsize=(8, 6))
sns.heatmap(confusion, annot=True, fmt='d', cmap='Blues', cbar=False)

# Add labels and title
plt.xlabel('Predicted Labels')
plt.ylabel('True Labels')
plt.title('Confusion Matrix')

# Show the plot
plt.show()


from skimage import io, color, filters, feature

# Compute ordered gradients for each image
ordered_gradients = []
#original_image = grayscale_images.reshape(64, 64)
for img in grayscale_images:
    img = img.reshape(64,64)
    gradient_x = filters.sobel_h(img)
    gradient_y = filters.sobel_v(img)
    
    # Compute the gradient orientation (in radians)
    gradient_orientation = np.arctan2(gradient_y, gradient_x)
    
    ordered_gradients.append(gradient_orientation)


ordered_gradients[1]

array([[ 2.35619449, -1.64756822, -1.67046498, ...,  2.35619449,
        -1.03037683, -0.78539816],
       [ 2.67794504, -2.0647377 , -1.74546853, ...,  3.05093277,
         2.89661399,  0.        ],
       [ 2.9996956 , -2.2318395 , -1.46013911, ..., -2.94419709,
         2.46685171,  1.9513027 ],
       ...,
       [ 0.5404195 ,  0.78539816,  0.78539816, ..., -2.35619449,
         3.14159265,  2.35619449],
       [ 1.57079633,  1.37340077,  0.78539816, ..., -2.8198421 ,
        -2.8198421 ,  3.14159265],
       [ 1.89254688,  1.42889927,  0.78539816, ...,  0.        ,
        -1.89254688, -2.35619449]])


#Scaling the transformed data using min-max scalar to adjust any extreme data values in the dataset.

from sklearn.preprocessing import MinMaxScaler

transformed_ordered_gradients = []
for grad in ordered_gradients:
    temp = grad.flatten() # convert 2-D to 1-D
    transformed_ordered_gradients.append(temp)

scaler = MinMaxScaler()
transformed_ordered_gradients = np.array(transformed_ordered_gradients)

if len(transformed_ordered_gradients.shape) == 1:
    transformed_ordered_gradients = transformed_ordered_gradients.reshape(-1, 1)

# Fit the scaler on your data and transform it
scaled_images_grad = scaler.fit_transform(transformed_ordered_gradients)

# If you want to retain it as a numpy array (optional)
transformed_ordered_gradients = np.array(scaled_images_grad)


# Visualize the gradient orientations for some images
plt.figure(figsize=(12, 6))
for i, gradient_orientation in enumerate(ordered_gradients[:10]):  # Display the first 5 images
    plt.subplot(2, 5, i + 1)
    plt.imshow(gradient_orientation, cmap='hsv', vmin=-np.pi, vmax=np.pi)
    plt.title(f'Gradient Orientation {i + 1}')
    plt.axis('off')

plt.tight_layout()
plt.show()


#Apply KNN to verify the effeciecny of this technique
from sklearn.neighbors import KNeighborsClassifier
from sklearn.metrics import accuracy_score, classification_report, confusion_matrix
from sklearn.model_selection import train_test_split

#Use KNN to identify if the Ordered gradient method helps improve the clasification of vehicles and non-vehicles

# Split the dataset into a training set and a testing set (e.g., 80% train, 20% test)
X_train, X_test, y_train, y_test = train_test_split(transformed_ordered_gradients, labels, test_size=0.2, random_state=42)

# Create a K-nearest neighbor classifier
knn_classifier = KNeighborsClassifier(n_neighbors=5)  

# Train the classifier on the training data
knn_classifier.fit(X_train, y_train)

# Make predictions on the testing set
y_pred = knn_classifier.predict(X_test)

# Evaluate the classifier
accuracy = accuracy_score(y_test, y_pred)
confusion = confusion_matrix(y_test, y_pred)

# Print the results
print(f'Accuracy: {accuracy}')
print('Confusion Matrix:')
print(confusion)

Accuracy: 0.5351914414414415
Confusion Matrix:
[[1774    0]
 [1651  127]]


from skimage.filters import gabor_kernel
from scipy import ndimage as ndi
from scipy import stats

# prepare filter bank kernels

X = grayscale_images[:10]
Y = labels[:10]
kernels = []
for theta in range(8):# orientations
    theta = theta / 8. * np.pi
    for sigma in (1, 3, 5): # std
        for frequency in (0.05, 0.15, 0.25, 0.35): # frequency
            kernel = np.real(gabor_kernel(frequency, theta=theta,
                                          sigma_x=sigma, sigma_y=sigma))
            kernels.append(kernel)

            
# compute the filter abank and take statistics of image
# this can be drastically sped up using the 2D FFT
def compute_gabor(row, kernels, shape):
    feats = np.zeros((len(kernels), 4), dtype=np.double)
    for k, kernel in enumerate(kernels):
        filtered = ndi.convolve(row.reshape(shape), kernel, mode='wrap')
        _,_,feats[k,0],feats[k,1],feats[k,2],feats[k,3] = stats.describe(filtered.reshape(-1))
        # mean, var, skew, kurt
        
    return feats.reshape(-1)

idx_to_reconstruct = int(np.random.rand(1)*len(X))

gabr_feature = compute_gabor(X[idx_to_reconstruct], kernels, (h,w))
gabr_feature.shape

(384,)


%time 
# Apply the gabor filter on the selected sample data
gabor_stats = np.apply_along_axis(compute_gabor, 1, X, kernels, (h,w))
print(gabor_stats.shape)

CPU times: user 1e+03 ns, sys: 0 ns, total: 1e+03 ns
Wall time: 3.1 µs
(10, 384)


# Define the dimensions for each image
h = 6
w = 64

label_names = ["vehicle" if item == 1 else "non-vehicle" for item in Y]

# Iterate through each image vector, reshape, and display
for i, image_vector in enumerate(gabor_stats):
    # Reshape the image vector into the desired dimensions
    image_matrix = image_vector.reshape(h, w)

    # Display or process the individual image here
    plt.imshow(image_matrix, cmap='gray')  # You can specify the colormap
    plt.title(label_names[i])
    plt.axis('off')  # Turn off axis labels and ticks
    plt.show()


from scipy.spatial import distance

sample_images_gabor = np.array(gabor_stats)


# Calculate pairwise differences for sampled vehicle images
vehicle_distances = distance.pdist(sample_images_gabor, 'euclidean')
vehicle_heatmap = distance.squareform(vehicle_distances)


# Plotting
plt.figure(figsize=(12, 6))

plt.plot( 2, 1)
plt.imshow(vehicle_heatmap, cmap='hot_r', interpolation='none')  # Use 'hot_r' colormap for reversed colors
plt.title('Pairwise Differences for Vehicle Images')
plt.colorbar()


plt.tight_layout()
plt.show()

Lab Assignment Two: Exploring Image Data¶

Course: Maching Learning in Python¶

Semester: Fall 2023¶

Team: Perceptron¶

Authors¶

Data Sources:¶

Business Understanding:¶

Standard PCA¶

Randomized PCA¶

Standard PCA v/s Randomized PCA¶

Using Daisy for feature extraction¶

KNN to identify the effect of Daisy feature to the claissification model.¶

Exploring Gabor Filters and Ordered Gradient techniques:¶

Ordered Gradient¶

Gabor Filters¶

References¶